function [seed_coord] = seed_dclayer(nodiv,ptsperdiv,tho,thf,yo,yf,ae,be)
%   seed_dclayer  coordinates in dorsal column layer
%   seed_dclayer(nodiv,ptsperdiv,tho,thf,yo,yf,ae,be) splits the
%   dorsomedial aspect of the spinal cord into "nodiv" layers. This region
%   is known as the dorsal column (DC) of the spine.
%
%   Definitions:
%   nodiv = number of DC layers/divisions
%   ptsperdiv = number of locations per division
%   tho = angle (in degrees) of start of DC area
%   thf = angle ("") of end of DC area
%   [yo,yf] = y range spanned by DC area
%   ae = x axis of ellipse (ellipse in xy plane)
%   be = y axis of ellipse ("")

nopts=nodiv*ptsperdiv;
seed_coord=zeros(nopts,2); % 2 columns: x,y

% convert degrees to radians
tho=tho*pi/180;
thf=thf*pi/180;

% splitting arc azimuthally
thbnd=tho:(thf-tho)/nodiv:thf; % theta of each boundary between layers
thmid=( thbnd(1:end-1) + thbnd(2:end) )/2; % theta between layers
r_thmid=ae*be./sqrt((be*cos(thmid)).^2+(ae*sin(thmid)).^2); % elliptical rad.
x_thmid=r_thmid.*cos(thmid); % x coordinate
y_thmid=r_thmid.*sin(thmid); % y coordinate

% splitting down-up radius of ellipse
ybnd=yo:(yf-yo)/nodiv:yf;
ymid=(ybnd(1:end-1)+ybnd(2:end))/2; % midpoints

for ii=1:nodiv
    
    % displacement vectors
    ra=[0,ymid(ii)];
    rb=[x_thmid(ii),y_thmid(ii)];
    rc=rb-ra;

    for jj=1:ptsperdiv
        kk=ptsperdiv*(ii-1)+jj;
        f_jj=(2*jj-1)/2/ptsperdiv;
        seed_coord(kk,:)=ra+rc*f_jj;
    end
    
end

%% DEBUGGING
%   uncomment to visualize

% s  = 0:pi/100:2*pi;
% r_ellipse = ae*be./sqrt((be*cos(s)).^2+(ae*sin(s)).^2);
% x_ellipse = r_ellipse.*cos(s);
% y_ellipse = r_ellipse.*sin(s);

% hold on;
% clr_spec = input('');
% plot(x_ellipse,y_ellipse,'k-');
% plot(seed_coord(:,1),seed_coord(:,2),clr_spec);
% hold off;
